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Abstract 

Lattice Boltzmann Models (LBM) and Phase Field Models (PFM) are two of the most widespread 
approaches for the numerical study of multicomponent fluid systems. Both methods have been 
successfully employed by several authors but, despite their popularity, still remains unclear how 
to properly compare them and how they perform on the same problem. Here we present a uni- 
fied framework for the direct (one-to-one) comparison of the multicomponent LBM against the 
PFM. We provide analytical guidelines on how to compare the Shan-Chen (SC) lattice Boltz- 
mann model for non-ideal multicomponent fluids with a corresponding free energy (FE) lattice 
Boltzmann model. Then, in order to properly compare the LBM vs. the PFM, we propose a new 
formulation for the free energy of the Cahn-Hilliard/Navier-Stokes equations. Finally, the LBM 
model is numerically compared with the corresponding phase field model solved by means of a 
pseudo-spectral algorithm. This work constitute a first attempt to set the basis for a quantitative 
comparison between different algorithms for multicomponent fluids. We limit our scope to the 
few of the most common variants of the two most widespread methodologies, namely the lattice 
Boltzmann model (SC and FE variants) and the Phase Field Model. 

Keywords: phase field model, lattice Boltzmann, Navier-Stokes, Cahn-Hilliard, comparison, 
drop, leakage, spurious currents 



1. Introduction 

Since the beginning of the last century several approaches have been developed for theo- 
retical analysis and numerical simulation of multicomponent and multiphase flows. The most 
common theoretical approach to study such complex systems is based on the sharp-interface 
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assumption, in which the interface between the different fluids is considered of zero-thickness. 
Each component is characterized by its own physical properties, such as density, concentration, 
viscosity, etc. which are uniform over the portion of the domain occupied by the single fluid 
component. These physical properties may change in a discontinuous way across the interface 
and their jumps are determined by the equilibrium conditions at the interface. The evolution of 
the system is described by a set of conservation laws (mass, momentum, energy, etc.) separately 
written for each component. Their solution requires also a set of boundary conditions at the 
interface, which ultimately has to be tracked. The specific numerical difficulties involved with 
the interface tracking can be circumvented by the adoption of diffuse interface methods, like 
the ones focus of the present study: the multicomponent Lattice Boltzmann method (LBM) and 
the Phase Field Model (PFM). The drawback of these models is the high numerical resolution 
necessary to model real interface thickness that require a fictitious enlargement of the interface 
thickness. Despite this limitation, in recent years several researchers worked on the development 
and refinement of PFM as well as different variants of the multicomponent LBM. 

1.1. Introduction to LBM 

The kinetic theory for multicomponent fluids and gas mixtures has received a lot of atten- 
tion in literature UHZl]. Many of the kinetic models developed for the study of mixtures are 
based on the linearized Boltzmann equations, especially the single-relaxation-time model due to 
Bhatnagar, Gross, and Krook 18J, also named BGK-model: 

dfix, V, f) 

+ V ■ Vf{x, v,t) + a- Vyf(x, V, t) = Q(x, f) = (1) 

ot 

= -^[f(x,v,t)-f'^''^{x,v,t)], 

where f{x, v, t) is the probability density function to find at the space-time location (jc, t) a par- 
ticle with velocity v. The coUisional kernel, on the right hand side of Equation ([T]i, stands for 
the relaxation (with a characteristic relaxation time r) towards the local equilibrium f^''''\x, v, t) 
which, in turn, depends on the local coarse grained variables, as density and momentum: 



pix,t)^ Jf(x,v,t)dv pu(x,t)^ j f(x,v,t)vdv. (2) 

a ■ Vff(x, V, f) represents the effect of a volume/body force density, a, on the kinetic dynamics. 
Modern discrete-velocity counterparts of ([TJ, the so-called Lattice Boltzmann methods (LBM), 
are able to simulate multiphase and multicomponent fluids and have attracted considerable atten- 
tion from the scientific community 18]. The LBM is an discrete form of Boltzmann kinetic 
equation describing the dynamics of a fictitious ensemble of particles ] 19-22], whose motion and 
interactions are confined to a regular space-time lattice. This approach consists in the following 
evolution: 

Mx + c,Af, f + Af) - fi(x, t) = --[fi{x, t) - jf^ix, f)], (3) 

T 

where fi{x, t) is the probability density function of finding a particle at site x and time f, moving 
in the direction of the i-th lattice speed c, with / - 0, . . .,b. Systematic ways to derive the 
discrete set of velocities in these models are either the discretization of the Boltzmann equation 



on the roots of Hermite polynomials Il23l - |29[] or the construction of high-order lattices for more 



stable LBM based on entropic approaches H3(X 13111 . At the same time, the translation of the 



body/volume force a ■ Vv/(x, v, f) onto the discrete-lattice framework represented one of the 
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most challenging issues in the last years of Lattice Boltzmann resear chlllHlJ, 132-36.1 . Through 
one of the first approaches proposed in the literature, the so called Shan-Chen (SC) approach 



111 IL |12[1 . the non-ideal interactions have been introduced directly at the discrete lattice level 
among the constituent (kinetic) particles El [Hill. These lattice forces embed the essential 
features and are able to produce phase separation (i.e. a non-ideal equation of state and a non- 
zero surface tension) as well as a detailed diffuse interface structure. The application of the 
SC models has been particularly fruitful for many applications 1 3^ 3^ Nevertheless, its 



theoretical foundations have been object of debate in the recent vears i32[[37 , 40l, 41 1. mainly 
because of the thermodynamic consistency of the mesoscopic interactions involved. On the other 
hand, in the so called free-energy (FE) models 1,15-171, the collisional properties of the model 
have been chosen in such a way that the large scale equilibrium is consistent with an underlying 
free energy functional, embedding both hard core effects and weak interacting tails. In this case, 
more traceable theoretical foundations have been provided, at least from the point of view of a 



continuum theory 14311 . Among others, some studies have also performed where more elaborated 



lattice models, including the effect of an exclusion volume based on Enskog theory 132-341. 
effective equilibria [39J, or even effective SC forces, were designed to match the desired bulk 
pressure of a given fluid ll35i[36ll . 

1.2. Introduction to PFM 

The phase field model is based on the idea that the interface between two fluids is a layer 
of finite thickness rather than a sharp discontinuity. Across the interfacial layer the physical 
properties of the mixture components vary in a smooth and continuous way, from one fluid to the 
other This approach is based on the pioneering work of van der Waals |44], who first determined 
the interface thickness of a critical liquid-vapor mixture. In the PFM the state of the system is 
described, at any time, by an order parameter 4> - <l>{x), which is a function of the position vector 
X. The order parameter directly represents a physical properties of the fluid, such as its density, 
molar concentration, etc.; all the remaining properties are in turn modeled as proportional to (f)(x) 
ll45il43l . According to the diffuse-interface approach, the order parameter is continuous over the 
entire domain and it shows smooth variations across the interface between single fluid regions, 
where <p{x) assumes approximately uniform values. Coupling the continuous and diffuse repre- 
sentation of the system with a transport equation of the order parameter, the system evolution 
can be resolved in time. One of the best-known PFM is the Cahn-Hilliard equation 



an extension to binary mixtures of the work of Van der Waals ||44)1 . This equation is a transport 
equation for the order parameter, where the evolution of (pix) is proportional to the gradient of the 
chemical potential, fi. The chemical potential is defined in terms of the free energy functional, 

where the free energy, 'F[0], assumes suitable definitions according to the problem under anal- 
ysis (and also depending on which physical quantity has to be described) and is a conservative, 
thermodynamically consistent functional. The most common free energy formulation is given 
by the sum of an ideal part, 'FidW, and a non-local part, ;7^,/[V0]. The ideal part accounts for the 
tendency of the system to separate in two pure components and is derived from the thermody- 
namics of mixtures. The non-local part accounts for the diffusive interfacial region. Through the 
Cahn-Hilliard equation, the evolution of the order parameter is thermodynamically consistent and 



subject to a phase field conservation. As a result, the prediction of the interfacial layer does not 
deteriorate. In the case of density-matched fluid systems, the convective Cahn-Hilliard equation 
is coupled with a modified Navier-Stokes equation, where a surface tension (or capillary) forcing 
term, which is derived from the Korteweg stress |49], is introduced. This contribution yields 
to the Cahn-Hilliard/Navier-Stokes coupled equations system 150.- 321, w hich is also known as 
Model-H, according to the classification of Hohenberg and Halperin 05 3n . who studied the con- 
vective phase separation of a partially miscible fluid mixture. This model, originally developed 



to study critical phenomena L54.-,5 8. 
flows of Newtonian fluids 11461 D 



, was subsequently used by many authors to study two-phase 
59146111 . In this case, even if the fluids are in fact immisci- 
ble, molecular diffusion between the two species is allowed in the interfacial region. Thus, the 
thinner is the interface, the more realistic is the numerical solution. Realistic interface thickness 
requires high numerical resolution, which is usually beyond current computational limits. For 
this reason the interface is often kept larger than corresponding physical value (this approxima- 
tion holds also for lattice Boltzmann approaches). However, in spite of this approximation, the 
method has shown capabilities to capture complex interfacial dynamics in a wide range of real 



metnoa nas snown capaDUities to c; 
physical problems I45l 50l 62. 631. 



1.3. Advantages and disadvantages of the methods 

Multicomponent LBM and PFM have demonstrated excellent performances to predict the 
dynamics of multiphase and multicomponent flows. Yet, both methods show their own peculiar 
characteristics and drawbacks which can limit their use, performances and range of validity. A 
particular unexpected, and unwanted, feature of multiphase and multicomponent solvers is the 
manifestation of non physical velocities near equilibrium interface, present even for systems at 
rest. From a physical viewpoint the velocity should clearly vanish at equilibrium but, as it has 
been observed by many authors, small spurious currents most often exist in the proximity of the 
interfaces. In an attempt to remove these unwanted features several improvement to the LBM 
have been proposed in recent years ll64l - l69ll . It worths noticing that some of these improvement 
are capable to remove these spurious current to machine precision [67]. Spurious currents have 
also been observed in other numerical methods including the PFM [70-72J. 

Because of the magnitude of these spurious current drastically depend on the actual variant 
of the LBM or PFM, the comparison between the two methods may be somewhat ill defined. 
Here we aim at comparing the LBM vs. the PFM for their most widespread and used variants. 
Our answer will thus not provide a general statement valid for the two methods as such, but will 
still provide some extremely useful insight in what can be expected from the most employed 
variations of the approaches. As a side results, we will also quantitatively compare two of the 
most widely used lattice Boltzmann variants, the SC-LBM and the FE-LBM, under the same 
conditions (i.e. same diffuse interface model, same surface tension, same chemical potential, 
etc.). In order to achieve our goal, the problem of a one-to-one matching of the PFM with SC/FE 
multicomponent LBM needs to be addressed first. The one-to-one matching of the two methods 
gives also the opportunity to clarify how they compare with respect to the computational costs. 
In order to address these issues we start by analyzing the SC model for two population with 
inter-particle repulsion; the large scale continuum limit is reviewed and formulated in terms of 
a diffuse interface model with an underlying thermodynamic FE functional. In this way one 
of the crucial issues in the matching of SC model vs. the corresponding FE model is being 
solved. Then, starting from the matched SC/FE multicomponent LBM, a new formulation for 
the free-energy of the PFM is derived in order to directly compare them. Finally, a comparison 
of the numerical results obtained, on the same problem, from both LBM and PFM is presented, 
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focusing in particular on unwanted spurious currents or mass leakage in sheared suspensions. 
This work focuses only the case of binary mixtures, even if modeling more than two components 
is nowadays far from trial extension 173-771]. 



2. The lattice Boltzmann models 



2.1. The multicomponent Shan-Chen model 



In this section the multicomponent model introduced by Shan & Chen 111 U 11211 is reviewed. 
First the main properties of the model are recalled, then the equilibrium features (diffusive current 
and pressure tensor) relevant on the Jiydrodynamic scales are analyzed. Starting from a kinetic 
lattice Boltzmann equation 2^, 22 1 for a multicomponent fluid with A^, species iflslfTil] . the 



evolution equations over a characteristic time lapse Af read as follows: 

Ux + a At, t + At)- Mx, f) = - — IMx, t) - f^''\p„ u + T,F, /p,)] , (5) 

where fi^ix, t) is the probability density function of finding a particle of species s - 1, . . . , A^j. 
For the sake of simplicity, the characteristic time interval Af is assumed to be unity in what 
follows. The left hand-side of (|5]i accounts for the molecular free-streaming, whereas the right- 
hand side represents the time relaxation (due to collisions) towards local Maxwellian equilibrium 
/■'^'^*(p.s, m) on a time scale Tj JH [l^ 113. Elt] . The local Maxwellian is truncated at second order, 
leading to a sufficiently accurate approximation to correctly recover the hydrodynamic balance 
of an isothermal regime [, 19- 22.1. It should be noted that the equilibrium for the s species is a 
function of their local densities and common velocities, defined as: 

Mx,ty, U(x,t)^ ' ; — . (6) 



This common velocity receives a suitable shift from the force, Fs, acting on the s species 111 III 1311 . 
Fs may be either an external force or an internal force representing an intermolecular interac- 
tions: 

FAX, t) = -G,,,p,(x, 2 Z "/^^'(^ + f)'^'-' (7) 

where w|''*^ are suitable weights used to enforce isotropy of the resulting hydrodynamics j 2^ 
221. The weights, wf''\ are normalized as follows: 



^wf^^c?c5 = 5,,c^, (8) 

2 w'^^'^^clc'^C^c] = 4 {5atSsk + dakdbs + SasSbk) , (9) 

where 6ah is the Kronecker delta and - 1 /3 a constant. In the long-wavelength limit (i.e. 
when the fluctuations with respect to the equilibrium distribution function are small) the set of 
macroscopic equations associated with the kinetic model consist of the continuity equations (one 
for each component) and of the equation of motion for the total fluid momentum. Dealing with 
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the two species (i.e. A and B) under the assumption of the same characteristic time scale for all 
the components = tQ these equations are approximated by: 



dt 



(s) 



+ V ■ (p,h) = V ■ / 



-V ■ + V ■ (ri(T)Vu + ri(T)Vu' ), 



(10) 



(11) 



where p - 2.sP.s is the total density, u = Y^sPsUsIp is the baricentric (total) fluid velocity and 
P the pressure tensor with the property -V ■ P + Y(c^p) - YjsP- The kinematic viscosity is 
v(r) = (t - 1/2) and the dynamic viscosity is 77(7) = pv(r). The diffiisive current, J 
given by the following: 

J(A) ^ 



PaPb 
P 

_J(B) 



2/1 Pa Pb j \ Pa Pb j 



(12) 



and, for purely repulsive fluids, F^^'^^/pa^b yields: 



= -gABCs^PBA —CsVApBA + 

Pa,b ^ 



(13) 



where the coupling constant G, 
obtaining: 



gAB > 0. Both continuity equations may be subtracted. 



+ V ■ (0m) = V ■ 



(14) 



where = Pa - Pb and where 7*'^' = 2/*'**. The pressure tensor in the (global) momentum 
equation is: 



P = 



/ X 2 Sab . 2 Sab . 

(Pa + Pb) + SabPaPb + —PaI^Pb + —PbI^Pa 



(15) 



+ (VpA ■ Vpfi) 8 - clgAB^pA^PB + 



(T) 



where the extra spurious r-dependent contribution S'^'"' is small for t close to 1 /2. Its exact 
expression is reported in I42l] and, for the sake of simplicity, it will be neglected by setting t 
close to 1 /2, thus AT*'"' 0. Summarizing, the following equations system is found: 



1PaPbc\ I 1 

— r2 



p 

IpAPBT 
P 



YPa 

Pa 



YPb 

Pb 



(16) 



2 §AB 2 

SabCs ^(Pb -Pa)+ — VA(pb - Pa) 



, \ 2 Sab . 2 , 

(Pa + Pb) + SabPaPb + —PaI^Pb + —pb^Pa 

+ (VpA ■ Vpfi) d - CggAB^pA'^PB- 



(17) 



' Whenever timescales are different, the characteristic time maps directly onto an effective relaxation time. For the two 
species system (A,B) this assumes the form f = EAlAlEMlM.- similar expressions need to be used in the total baricentric 
velocity. 



The variables p - Pa + Pb and 4> — Pa - Pb can be introduced such as = (p + <;*)/2 and 
Pb — ip - <P)/'2- Under the assumption of an incompressible fluid at constant p, the derivative 
with respect to the density can be neglected in the pressure tensor: 



48AB 



P+—(P -0 — |V0| 



while the diffiisive current is mapped into the following expression: 



T 

2 



Vlog 



(T) 



or, in a more compact form: 



= M(0)V//(0) 



^Vlog 

2 ^ 



2 



Into the equation above, the density-dependent mobility has been introduced: 



and the chemical potential yields: 

P + 



P-, 



^-2 ' 



where the modified (through the relaxation time t) coupling constant reads: 



(18) 



(19) 



(20) 



(21) 



(22) 



(T-i) 



gAB, 



(23) 



which reduces to gAB in the limit t » 1/2. Summarizing, the set of incompressible equations 
which are approximated on a large scale are given by 



^ + V ■ (0H) = V ■ /(' 
at 



p -+«.V« 



-V ■ P + V ■ (ri(T)Vu + ri(T)Vu' ), 



(24) 



(25) 



where P and 7*'^' are given in ( l20b and ( fTSl ). It worths noticing that, for the models here consid- 
ered, the phase diffusion is dictated by the same relaxation time ruling the kinematic viscosity, 
i.e. fixed Schmidt number is considered. In order to decouple these two time scales a different 
more elaborated formulation based on multiple relaxation time schemes (MRT) 1 78 , 7^ should 
be adopted. 
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2.2. The multicomponent free energy model 

The other lattice Boltzmann variant that we consider in this work is the Free Energy (FE) 
based model lll5l - ll7ll . Within this approach the equilibrium properties of the model are con- 
structed to be consistent with an underlying continuum free energy functional of the order pa- 
rameter, whose formulation is the following: 



d^:, (26) 



where V{4>) is the bulk free energy density and k a constant parameter associated to the surface 
tension at the interface. The chemical potential /i(0) and the pressure tensor are derived from 
thermodynamics identities: 

... 6TW dVicp) 



P 



(^PbW - Kcf>A(l> - ^1 V0|-) 6 + kV(I>V(I>, (28) 

where Pb(<p) - (pdV{(p)/d(p - V((p) is the bulk pressure. For the out of equilibrium properties, the 
mesoscopic dynamics is given by the following two coupled kinetic equations: 

fi(x + cAt,t + At)-Mx,t) = --mx,t)-fl"%,<p,u)], (29) 

T 

gi(x + cAt,t + At)-gi(x,t) = - — [gi(x,t)-gf''\4>,u)], (30) 
where the macroscopic density and velocity are given by: 

p(x,t)^yMx,ty, u(xj)^^^^^p^. (31) 

The order parameter density, (p, is obtained from the populations, gc. 

(l>ix,t)^Yjgiix,t). (32) 

The equilibrium distribution function is then chosen with a polynomial form in the kinetic veloc- 
ities, so that: 

Y,f'^CiCi^P+puu, (33) 

Y^gf''>CiCi^SriJ + 4>uu, (34) 

where T controls the order parameter mobility M = F^t^ - 1/2^0. The previous kinetic model 
approximates the following macroscopic equations in the incompressible limit: 



^ + V ■ (0H) = MAfi{<f>), (35) 
ot 



du 



p\^^ +uVuj^ -V ■ P + V ■ (77(t)Vh + 77(t)Vh^), (36) 

where the viscosity is the same as in ( l25T l. The presence of extra spurious terms in the large scale 
expansions (|35])-(|36]| has been discussed in fSOl]. Through the numerical simulations presented 
in this work, these terms are directly evaluated and found to be negligibly small. 



in the free energy approach the mobihty can be easily changed by properly tuning both Tg and T. 
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2.3. Matching in the long-wavelength limit 

In this section the details of the matching between equations (l24l i-(l25Tl and (|35]|-(|36]| are 
reported, representing the large scale limits of the SC and FE models, respectively, when com- 
pressible effects are negligible (p constant). Since the transport properties in the hydrodynamical 
equations are in fact the same, the crucial issue for the matching is to set the same thermody- 
namical background functional. First the continuity equation (l24l l is rewritten in a formulation 
similar to equation (l35l) : 



^ + V ■ (0h) = V ■ j'-f^ = V [M(0)Vju] , 
ot 

where the chemical potential can be derived from the following r-dependent functional: 



(37) 



(38) 



(T) 2(P 

Saba's — 



(t) 4 



from which, by definition, the chemical potential is extracted: 



- V- 



drw 



= flog, 
2 \p 



p + i 



W) 

°AB 4 A JL 



(39) 



and it coincides with equation ( l22l l. In order to set an accurate matching between the two lattice 
Boltzmann approaches, two more issues must be analyzed. Although one may think to integrate 
the FE model (|29Tl-(l30ll with the functional ( l38T l instead of (|26] |. both continuity equations are not 
exactly the same due to the (fi dependency of the mobility in equation ( l37b . However, equation 
(IJTI i can be approximated with the mobility in the center of the interface (that means negligibly 
small density changes <k p): 



(40) 



M(0,p,T)*M(p,r)=p|T-- . 



With this approximation, the advection-diffusion equation for the order parameter (p obtained 
from the SC model and from equations ( l30b and ( l34l i of the FE model with the functional (l38T l 
are set to be the same in the hydrodynamical limit. Furthermore the pressure tensor should be 
matched between the two (FE and SC) LBM descriptions. The functional of equation ( l38l l con- 
tains a spurious and nonphysical r dependence hidden inside the definition of g^'^ = (jgAB)/(j - 
1/2). Only towards the limit of large t the functional ( l38l l becomes independent of t, whereas, 
for finite t, this extra spurious dependency should be kept inside equation ( |38] ) in order to ac- 
curately match both SC and FE descriptions. This difference emerges in the SC model because 



the collisional kernel does not conserve momentum locally 11131 11411 . Therefore the lattice mo- 
mentum should be shifted by a factor FAt/2 to reproduce the hydrodynamical equations. It is 



^Equation 1371 has been simulated with and without the constant mobility by adding proper counter-terms to the 
lattice Boltzmann equation. No substantial differences have been observed for the features analyzed in this paper. 
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important to remark here that such a spurious contribution may be removed by curing the discrete 
forcing effects as described in 181,] . This technique has already been adopted in recent simula- 
tions with the multicomponent SC model in the framework of the so-called multiple relaxation 
time schemes (MRT) IITSl pTglPl Differently with respect to the chemical potential, the pressure 
tensor of the SC description can be directly extracted: 



2 l P + <l> \, ( P + <P \ . 2/P-0\, ( P-^ \ 



(41) 



gABcj 



that is the large t limit of (l38T l. In fact, following standard procedures, the conserved current 
(pressure tensor) under the hypothesis of translational invariance of the free energy can be de- 
rived. The general expression is: 



(42) 



where the mass conserving free energy has been used: 



r'"n'i>,p, V0] = r[c/>,p, v,p]-Ai4>- a2p. 



(43) 



Ai^2 are two Lagrange multipliers introduced to ensure the global conservation of (p and p. Their 
values are readily evaluated using the Euler-Lagrange equations (should be noticed that there are 
no gradients of p in the free energy functional): 



d(p 



this yields: 



^1 = 



dV{(/>,p) 8ab4 



dp 



A(f), 



d(p 4 

At the end, the following pressure tensor is obtained: 



^2 = 



5(Vp) 

dV{<p,p) 

dp ■ 



0, 



p = 

+ 



, dV{4>,p) dV{<f>,p) 

a P' 



dp 



-V(cf>,p)-K<pAcp- -\Ycl>\' 



where k - (c'^gAB)/^) and the bulk potential V(<p) given by: 



gABCs ■ 



(44) 



(45) 



(46) 



(47) 



' Since the goal of the present work is to propose guidelines to match the SC and the FE approaches we limited the 
analysis to the basic versions of the algorithms. 
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It can be immediately noticed: 



^ — 0^ — +P — y('P,P)^c^p + Cs—(p-<p), (48) 

and, at the end: 

P = (cIp + 4 ^(p2 - ^2) - ^0A,^ - ^|V,^|2) ,5 + KV<pVcf>, (49) 

that coincides with ( fTSl l. We conclude that the momentum equation from the SC model and from 
equations ( |29] l and (l33l l of the FE model with the functional i41i are set to be the same in the 
hydrodynamical limit. 



3. The phase field model 

In this section the phase field model for a multicomponent fluid system is analyzed. Starting 
from the features of the PFM proposed by Badalassi et al. [46] and Yue et al. [61J, a new 
formulation of the free energy functional is derived to match the SC/FE multicomponent LBM 
reported in Section|2] For an isothermal binary mixture, where is the relative concentration of 
the mixture components, the free energy functional reads: 



dx, (50) 



where V(^) is the ideal part of the specific free energy, (K/2)\Vcf>\^ is the non-local part of the 
specific free energy and a: is a constant positive parameter associated to the surface tension at the 
interface. Recalling equation (|4]l the chemical potential fi(<p) is: 



64> d(j) 



dV{(h) 

-^-kA(P. (51) 



For unsteady problems, the evolution of the system is described by the Cahn-Hilliard equation 
[47, 48.1 , in which the evolution of (p(x) in time is proportional to the gradient of the chemical 
potential: 

^ + H ■ V0 = V ■ [M(^, p, T)Vfi(<p)] . (52) 
ot 

According to Section l23l the functional chosen to describe the free energy of the system, is the 
T-dependent 'F^'^\(p] of equation ( l38l l. Thus the chemical potential ju*'^'(0) reads: 



c: 



where the surface tension parameter is a- = Sab'^s ^^'^ '■^^ coupling coeflicient is g^^^ - 
{tSab) / (t - 1 /2). Substituting equation ( l53T l for the chemical potential of equation ( |52] | and con- 
sidering the mobility parameter M uniform over the domain, the Cahn-Hilliard equation yields: 

— +u-Vcf)^ M(p, T)A^'^\(f>). (54) 
ot 
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Focusing on density-matched (p ^ const.), viscosity-matched (y ^ const.) and incompressible 
(V • H = 0) Newtonian fluids, the flow field evolution is described by the modified Navier-Stokes 
equation 14611 : 

— + H ■ Vh = -Vp + v(r)V ■ (Vh + Vh^) + F^, (55) 
at 



where is the local capillary stress due to the concentration field 118211 . In order to match the 
same hydrodynamical properties of the multicomponent LBM, F^ is modeled according to the 
results of Sectionl23] 



F^ = ju(0)V0, (56) 

where ju(0) reads: 

Equation ( fSTT i is obtained applying the definition ^ to the functional (HTt . which is derived from 
the matching procedure of Section l273] Within this formulation, /i(0) represents the large t limit 
of the equation (|53] ). and //^'"^ — > ju when r » 1/2. Finally, the PFM equations implemented in 
this work are the following: 

V ■ M = 0, (58) 

^ + H ■ V(/. = M(p, t)Aii'^\(P), (59) 
ot 

du r 

— +u Vu = -Vp + v(t)V ■ (Vm + Vu' ) + fi((f>)'V(f), (60) 
ot 

where the chemical potentials p and are given by equations ( |57] | and (|53|) respectivelv. Their 
formulation is similar to that used by Mauri et al. [54], Vladimirova et al. 05511 and Molin et al. 
ll56 1 . The mobility coeflicient of equation ( |59] | is M(p, r) = p (t - 1 /2), whereas the kinematic 
viscosity of equation ( l60l l reads v(r) = c^(r - 1 /2). The coefficients t, p, gAs and = 1/3 are 
the LBM input parameters and their definitions are reported in Tabm Equations (I58]l-(l60ll have 
been rewritten in a dimensionless form using the following dimensionless variables: 

u X t p <p 

u = — , jc" = — , r = , p' = — — , 0" = — . (61) 

U H H/U ^ pU^/H ^ (j)* ' 

where H, U, (p* are the characteristic length, velocity and concentration, respectively. Starting 
from equation (|53] |. the chemical potential yields to: 



1 ,„JP" + 

( 

consistently equation (ISTT i reads: 



7;^logl^— ^l-r-C/i^A^-, (62) 
Gab \p 



= -\ log I -r- Ch^Acf,-, (63) 

where p" - p/(p*. Gab - gAB<P* and G^'^ = g^AB'f'*- "^^^ dimensionless form of equations 
(|59]i and ^ are: 

V ■ = 0, (64) 
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Yt 



— — +u ■ V0 



du 



— — + H ■ Vm = -Vp 



Re 



Pe ^ 
(Vu + Vh"^) + 



1 



(65) 



(66) 



ReChCa 

The non-dimensional groups introduced into the equations system (l62ll-(l66ll are the Cahn num- 
ber, the Peclet number, the Reynolds number and the capillary number, which are defined as 
follows: 

l_ _ ^ vpU 

H' 



Ch 



Pe 



Re 



Ca 



(67) 



The Cahn number is the ratio between the interface thickness ^ and the length-scale H, the Peclet 
number is the ratio between the diffusive time-scale H^/(J3^'^^M) and the convective time-scale 
H/U. The Reynolds number is the ratio between the inertial forces HU and the viscous forces v. 
Finally the capillary number is the ratio between the viscous forces vpU and the capillary forces 
at the interface. All the dimensionless groups and the definitions of the fluid properties with 
respect to the LBM input parameters are reported in TabQ] The equations system (l62b-(l66t can 
be easily extended to systems with different viscosity and different mobility between the two 
components 146 1 . Mismatched-density systems require the adoption of a density -based PFM 



3.1. The numerical method 

The equations system (|64]|-(|66]| is directly solved using a pseudo-spectral algorithm in which 
the field variables are transformed into wave-number space by means of Fast Fourier Transform 
(FFT). The non-linear terms are evaluated in the physical space an then re-transformed into the 
wave-number space (Soldati and Banerjee 1831]). In order to reduce the time-step constraint 
required by the Cahn-Hilliard equation, an operator splitting technique has been introduced on 
to equation ( |65] |. The procedure adopted is similar to that discussed in [46.1 and it yields (where 
the apex "-" has been removed for sake of simplicity): 



dt 



-u -VS + —A 

^ Pe 



1 



log 



p-< 



(68) 



— iCh' 
Pe 



\a^<P+ — 

2Pe / Pe 



2A^ 



At 
iFe 



By explicit treatment of the first term on the right-hand-side of ( l68l l and implicit treatment of the 
second term on the right-hand-side of equation ( |68] |, an efficient semi-implicit discretization is 
obtained. The time step advancement for both equations (|66] | and ( l68T l is obtained using a two- 
level Adams-Bashfort scheme for the explicit non-linear terms, and the Crank-Nicholson scheme 
for the imphcit linear terms. Applied to equation ( |68] | this implicit-explicit combination reads: 



(69) 



At 2 

where the explicit S term and the implicit *F term are: 



-u ■Vd>+ —A 

^ Pe 



Pe \ 2Pe I 
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1 

2 J. 



log 



P-' 



3(p 



(70) 



More details on the numerical algorithm can be found in ll83ll . 



(71) 



4. Numerical tests 

In this section the numerical results obtained from both the PFM and the LBM approaches 
are discussed. Two different tests have been performed with the three models, SC-LBM, FE- 
LBM and PFM, as described in Section ITTl Section l272l and Section [3] respectively. First, the 
equilibration of a two dimensional static droplet has been simulated until the steady state has 
been reached. Then, starting from the settled droplet, its deformation under a Kolmogorov flow 
has investigate. For the numerical analysis a SC-LBM model with r = 0.55, p = 1.4 and 
gAB = 0.164 has been used. The interaction parameter r has been chosen small enough to avoid 
spurious contributions from the pressure tensor. The other parameters (gAB and p) have been 
chosen in order to obtain (fi - pA - Pb - ±^ inside the pure components. The simulations have 
been carried out on a two dimensional fully periodic grid of 100 x 100 nodes, on a computational 
domain of dimensions x L,, - 100 x 100. The same simulations performed with the SC model 
have been repeated with the FE-LBM. Then the PFM simulations have been performed on a two 
dimensional fully periodic gricfl composed by 128 x 128 nodes on a computational domain of 
dimensions Lj x L,, = 2nH x 2nH - 100 x 100. In particular, PFM dimensionless numbers 
have been calculated following the definitions reported in Tab [Tl and their values are collected in 
Tab [2] In order to match = ±1 in the regions of pure components, the scaling concentration 
parameter has been chosen 0* = 1 . As result SC/FE-LBM and the phase field model have been 
set with the same interface and transport properties and thus the results can be both qualitatively 
and quantitatively compared. 

4.1. A steady droplet 

A two dimensional static droplet with radius approximately /?/Lv - 1/4 initiated in a rest- 
ing fluid has been studied. The simulations have been run letting the drop attain a stationary, 
equilibrium state. The kinetic energy at the curved interface (Fig[T]l and the associated station- 
ary configuration of velocity field and order parameter (Fig [2] and FigO have been measured. 
The spurious currents shown by the PFM simulations are found almost one order of magnitude 
smaller that the ones from the LBM. On the other hand both EE and SC LBM models show 
comparable velocity magnitudes. The flow fields patterns are similar within the three methods, 
even if in the PFM case these currents are confined in a thinner layer along the interface. Fig [4] 
displays the temporal evolution of the surface tension cr, which is proportional to the difference 
between the pressure inside the drop (P,) and the pressure outside the drop (PoY- 

<T{t) = R [P,(f) - Poit)] ■ (72) 

Both the left panel of Fig[3]and Fig[4]show that the interface properties of the three models are 
the same (i.e. the same interface structure and the same surface tension at the curved interface). 
Nevertheless, the oscillations observed in the LBM models in Fig [4] are probably nonphysical 
pressure fluctuations which are ruled out in the PFM simulations. The most important feature 



''Due to the use of a spectral solver it is convenient to choose a power of 2 for the number of grid nodes. 
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observed appears to be the change in magnitude and structure of spurious currents: in both 
the LBM simulations approximately the same spurious currents are found while their intensity 
is reduced in the case of PFM simulations. Nevertheless, it is important to remark that the 
LBM used here are basic versions of the two widely adopted approaches. Improvements can be 
obtained by curing discretization errors in the computation of the intermolecular force causing 
parasitic currents as described by Lee & Fischer ll67ll . Such improvements have been shown 
to eliminate currents to roundoff if the potential form of the intermolecular force is used with 
compact isotropic discretization. For the sake of completeness the results obtained with the Lee- 
Fischer scheme have been reported in Fig [31 

4.2. Droplet deformation under Kolmogorov flow 

Starting from the equilibrium droplet obtained from the simulation of Section 14.11 a sinu- 
soidal forcing has been applied on the flow field until a new stationary state was reached. The 
forcing term has been chosen with the following formulation in order to generate a Kolmogorov 
flow: 



The evolution of the total kinetic energy reported in the right panel of Fig [7] shows a good 
matching between the three models, thus the same hydrodynamical transport properties have 
indeed been imposed. Little discrepancies (less than 10%) between PFM and LBM simulations 
are shown in the total kinetic energy of the steady state. Moreover little deformation differences 
can be observed for the concentration isocontours of Fig [6l with the PFM showing a slightly 
less deformed drop with respect to both LBM. On the contrary both LBM models (SC and FE) 
show the same concentration profile. The differences in kinetic energy and deformation seem to 
be consistent one with the other, in fact the less deformed the droplet, the less energy is absorbed 
from the flow and thus the higher the total kinetic energy (for same external forcing). Qualitative 
snapshots of the kinetic energy are reported in Fig[3] where the isocontours of = -0.9 and - 
0.9 have been superposed to show the interfacial layer location. Similar patterns and magnitudes 
are observed within the three models, confirming the correct matching of the models. To test the 
importance of spurious mass flux across the interface, the mass leakage L has been monitored in 
time: 



where M(f) and Mq are the number of computational nodes inside the droplet at time t and at 
time f = 0, respectively. The nodes with a local order parameter ^{x) > (pT, with (pj - 0.0, have 
been considered as belonging to the droplet. The results reported in Fig [8] show a negligible 
mass leakage (less than 0.5%) for the PFM and FE, whereas it slightly higher (in the order of 
l-3%)forSC-LBM. 




(73) 



where A = 10 In the PFM model the dimensionless forcing term is: 




(74) 
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5. Conclusions 

In this work two of the most widely adopted approaches for the numerical study of multi- 
component fluid systems, the Lattice Boltzmann Models (LBM) and Phase Field Models (PFM) 
have been compared on equal footing. First the Shan-Chen (SC) multicomponent LBM and the 
Free Energy (FE) LBM have been reviewed and analyzed. Focusing on the specific case of phase 
separating fluids with two species, the long-wavelength limit (i.e. the hydrodynamical limit) of 
both lattice Boltzmann models has been reviewed and a criteria to match the two models has 
been developed. Then, on the basis of the LBM matching, a new formulation for the free-energy 
involved into the PFM has been derived. In this way the analytical matching between PFM and 
LBM has been obtained. Finally the three models (SC-LBM, FE-LBM and PFM) have been 
numerically tested against controlled benchmarks with both steady and moving interfaces. Three 
main advantages emerge from this kind of analysis: first, from the theoretical point of view, it 
has been verified that the PFM equations are well approximated by the large scale limit of the 
multicomponent LBM. In fact, the convergence of the lattice Boltzmann models towards the 
diffuse-interface hydrodynamics may be questionable due to the presence of steep interfaces, 
where the local gradients are high enough to prevent the usual Chapman-Enskog analysis to 
be applied |20]. Second, from the computational viewpoint, some of the undesirable features 
emerging in the LBM simulations have been shown to exist also in the PFM, but these appear 
somewhat reduced, at least when basic versions of LBM are considered. Third, this study offers 
the possibility to test the performances of the different methods simulating the same physical 
system. The PFM shows a computational cost almost three times higher to the analogous LBM. 
On the other hand the quantitative results of a PFM appears to be more accurate. We believe 
this study helps clarify important issues beyond the choice of the either a lattice Boltzmann or a 
phase field based approaches for multicomponent fluid dynamics. 
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M = p(T-i) 



_ gABC^ 

~ 4 



v(t) 

Table . 1 : Definition of the Lattice Boltzmann input parameters and of the corresponding phase field model dimensionless 
groups. 



p 

1.4 


T 

0.5 


gAB 

0.164 


8ab 
1.804 


1 

3 


Pe 
0.7930 


Ch 
0.0256 


Ca 
0.0021 


Re 
1.000 





Table .2: Definition of the Lattice Boltzmann and phase field input values. 
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Figure .1: Contour plots of the local kinetic energy, per unit density, ^(m^ + uj) due to the spurious currents for the FE- 
LBM (top), SC-LBM (center) and PFIVI (bottom) methods when simulating a stationary two dimensional droplet. The 
snapshots are taken at the same time when a steady state has been attained. The intensity and structure of the spurious 
kinetic terms comparable in the lattice Boltzmann models and is reduced by a factor 100 in the PPM model (i.e. a factor 
10 on the velocity magnitudes). 
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Figure .2: Vector plots of the velocity field due to the spurious currents contributions for the FE-LBM (top), SC-LBM 
(center) and PFM (bottom) simulations for the two dimensional static droplet. The plots are taken at the same time when 
the steady state has been reached. The velocity field of the Lattice Boltzmann simulations (top and center plots) have 
been magnified by a factor 10"* whereas the vector field of the PFM simulation have been magnified by a factor 5 ■ Iff* 
for the sake of readability. 
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Figure .3: Local concentration and velocity profiles from FE-LBM, SC-LBM and PPM simulations of a two dimensional 
static droplet. On the left the local order parameter, if) = p^- Pb, is plotted as a function of the coordinate, v, for fixed 
,v = 50. On the right the local (spurious) velocity in the vertical direction, «,.. is plotted as a function of the coordinate, y, 
for fixed .v = 50. The order of magnitude of both spurious contributions is comparable in the lattice Boltzmann models 
while is reduced by a factor 10 (for the velocity, 100 for energy) in the PFM model. Improvements in the LBM can be 
obtained by curing discretization errors in the computation of the intermolecular force as described by Lee & Fischer 
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Figure .4: Time evolution of the surface tension (Laplace test) from FE-LBM, SC-LBM and PFM simulations of a two 
dimensional static droplet. Starting from the same initial conditions, the local value of surface tension tr(f) is plotted as 
a function of time. 
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Figure .5: Contour plot of the local kinetic energy per unit density ^(Hj + hJ) for the FE-LBM (top). SC-LBM (center) 
and PPM (bottom) simulations of a two dimensional droplet deformation under Kolmogorov flow. The plots are taken 
at the same time when the steady state has been reached. Similar magnitudes and patterns can be observed for all the 
models. 
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Figure .6: Isocontour plot of the concentration field at = for the FE-LBM, SC-LBM and PPM model for a stationary 
two dimensional droplet under shear. 
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Figure .7: Time evolution of the total kinetic energy from FE-LBM, SC-LBM and PFM simulations of a two dimensional 
droplet deformation under Kolmogorov flow. Starting from the same initial conditions, the total value of the kinetic 
energy //("^ + "y)'l^dy is plotted as a function of time. Similar evolution in time is registered for all the models, even if 
the PFM showed an asymptotic value higher than the LBM. The CPU elapsed time of both PFM (TpfM) and SC (Tsc) 
methods have been measured through this simulation. The Lattice Boltzmann method was roughly three times faster 
than the analogous Phase Field Model = 2.9). 
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Figure .8: Time evolution of the relative leakage of the order parameter from FE-LBM, SC-LBM and PPM simulations 
of a two dimensional droplet deformation under Kolmogorov flow. Starting from the same initial conditions, the relative 
leakage L( t) is plotted as a function of time. Longer simulations show that the leakage reached the saturation value for 
all the three methods. 
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